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Abstract 



When observing from the ground one has to deal with the turbulence from the 

atmosphere and the turbulence caused by the telescope's surroundings. In or- 
der to find the best sites for future solar telescopes and develop multi-conjugate 
adaptive optics (MCAO) for them, the contributions to seeing have to be char- 
acterized for heights up to at least 12 km above the telescope. 

This characterization is done with a method that measures differential image 
displacements using several subapertures. By extending the ordinary DIMM- 
and S-DIMM-methods to measure displacements at different field angles, the 
height distribution of seeing can be measured. This extended method is called 
S-DIMM-I-, where the letter S denotes the use of the Sun as source and the -|- 
signifies that it is an extension of the former methods. 

Observations are made with the Wide Field Wavefront Sensor (WFWFS) 
mounted at the Swedish 1-m Solar Telescope (SST) at La Palma. The re- 
sults from these measurements are compared with the Shadow Band Ranger 
(SHABAR), also mounted on the SST. 

The first results of these comparisons show good correlation between the two 
instruments for heights up to 500 m. Results from the WFWFS also show 
that the dominant seeing at this site comes from layers close to the ground and 
at high altitude and that Pried's parameter, ro, is larger than 40 cm for the 
intermediate layers. 
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st Chapter 



Introduction 



Image quality is very important for all types of observations and it has been 
known for quite long that the quality of an image observed with a telescope 
does not only depend on which type of telescope that is used, the site where the 
telescope is located matters. The atmosphere, through which observations are 
made (with ground based telescopes), perturbs the wavefronts causing image 
degradation. 

Different methods exist to compensate for the wavefront distortion caused by the 
atmosphere. The nature of the distortion must be known in order to compensate 
for it. Different ways of quantifying these distortions exists as well. 

The distortion is different for different sites, making the location of the telescope 
very important. The distortions are also changing rapidly with time. Methods 
for quantifying this distortion and for compensating it in real time are therefore 
needed. 

A brief introduction to the origin of these distortions and how to compensate 
for them are given below together with a short description of the next European 
solar telescope. 

This chapter ends with a short outline of this report. 



The atmosphere is not, as one might imagine, homogeneous, it is turbulent and 
consists of turbulent eddies with different temperatures. The index of refraction 
is primarily a function of temperature and it will therefore be different for the 
different layers in the atmosphere. Light is refracted when traveling through the 
turbulent atmosphere and the optical path length for one line-of-sight will differ 
from another. Plane wavefronts will be disturbed after passing through the 
atmosphere making the observed objects blurred, destroying the visible details. 

One example of this is stars that seem to twinkle when you observe them a clear 
night with your naked eye. The twinkling is not due to the stars themselves, 
it's caused by high-altitude turbulence in the atmosphere and the effect is called 
seeing. Figure [TTT] schematically describes this. 



Seeing 
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These rippled wavefronts result in different forms of image degradation. Light 
that follows different paths through the atmosphere can form separate images 
of the same object. Labeyrie [13] referred to this as speckle patterns. 




Figure 1.1: Plane wavefronts enter the atmosphere and are distorted by the turbulent 
layers. The different objects in the turbulent layer represent regions with different 
refraction indexes n. A sample light ray shows how refraction in these objects leads 
to the perturbed wavefront. 



The turbulence that causes bad seeing can be anywhere, from the air inside the 
telescope up to the high atmosphere. High-layer seeing, Fig. 1 1.2a] is significantly 
less intense than ground-layer seeing, Fig. |1.2b[ for daytime observations. The 
different layers will cause different image distortions. High-layer seeing causes 
differential distortions because different parts of the field of view (FOV) will 
cross different regions of air at high altitude. 





Aperture 

(a) 




(b) 



(c) 



Figure 1.2: (a) High-layer seeing. Different parts of the atmosphere affect different 
wavefronts so that the distortion varies across the field of view, (b) Ground-layer 
seeing. The field of view is affected equally since wavefronts from all directions traverse 
through the same turbulence, (c) Multi-layer seeing. It is most common that seeing 
occurs both in higher layers and closer to the ground. High-layer seeing is normally less 
intense than ground-layer seeing. Pictures are adapted from T.I.M. van Werkhoven. 
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1.2 Seeing correction 

Short exposure images need to be taken in order to observe speckles, otherwise 
they wiU be averaged out. Long exposures will cause overall image degradation 
since the speckle patterns will overlap. The resolution of the image will therefore 
be degraded. In 1970, Labeyrie [13] presented techniques to partially restore the 
original image in post processing from a burst of short images. A disadvantage 
with this method is that it is limited to good seeing and to small telescope 
apertures. 

Another method is adaptive optics, which is the most powerful technique for 
all astronomical applications [1]. This technique uses a wavefront sensor that 
controls an optical component, a mirror whose surface can be deformed, in- 
troducing a controllable counter-distortion to the wavefront. The sensor will 
measure the distortion of the incoming wavefront and feed this information to 
the mirror. The mirror will then correct the distortion in real-time and the large 
scale structure of the original wavefront will be restored. Figure |1.3| schemati- 
cally describes the corrector in this technique. 




Figure 1.3: The correcting device of an adaptive optics system, a deformable mirror. 
The deformable mirror is deformed so that the wavefront distortions caused by the 
atmosphere are compensated for. 



By the use of several deformable mirrors, a Multi-conjugate Adaptive Optics 
(MCAO) system can be built. The deformable mirrors conjugate different alti- 
tudes and the turbulence can be compensated for in a three-dimensional way. 

A more detailed description of adaptive optics is given in Chapter [3j 



1.3 European Solar Telescope 

The European Solar Telescope (EST) is a future 4- meter solar telescope that will 
be located in the Canary Islands. The project involves institutions and partners 
from 15 European countries. The project is divided into different workpackages 
and one of them (WP08000) is responsible for the site characterization [7]. The 
turbulence above the observatories of the Canary Islands are studied in order 
to achieve turbulence profiles for the MCAO system planned for the EST pO] , 

The site will be characterized with two different instruments, a long-base 



SHABAR (see Sec. 4.3) and a wide-field wavefront sensor (see Sec. 3.2.2), op- 
erating at the two sites, the Observatorio del Teide on Tenerife and the Obser- 
vatorio del Roque de los Muchachos on La Palma [7]. One SHABAR and one 
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WFWFS are placed on the Swedish 1-m Solar Telescope on La Palma and the 
other WFWFS and SHABAR will be placed on VTT (the Vacuum Tower Tele- 
scope) on Tenerife. The advantage of placing both a WFWFS and a SHABAR 
on the same building is that the seeing measured by the both instruments can 



be compared [21] • This comparison is done in section Sec. 6.3 



The aim of the observations and the data analysis is to understand the isopla- 
natic patch that is attainable with the MCAO and to compare La Palma's and 
Tenerife's high-altitude seeing [7]. 



This report will mainly cover the data reduction and analysis from the WFWFS 
mounted on the Swedish 1-meter Solar Telescope (SST) and comparison with 
the SHABAR, also mounted on the SST. Chapter [2] gives a theoretical back- 
ground to seeing and atmospheric turbulence. Chapter [3] describes the different 
components of an adaptive optics system. Chapter [4] describes different meth- 
ods to quantifying seeing and ends with the optical setup at the SST. Chapter [s] 
deals with the data processing made after data were taken. Chapter |6] provides 
some initial results and comparison between the WFWFS and the SHABAR 
mounted on the SST. In Chapter [7] some conclusions from these results are 
drawn and some future aspects are discussed. 



4 




nd Chapter 



Quantifying seeing 



The characteristics of the wavefront aberrations that need to be compensated for 
must be weU known in order to design a good AO system. These aberrations are 
random and can therefore only be described statistically. The statistics describe 
the seeing conditions, which evolve with time. In order to design MCAO (and 
other instruments) the statistics of their evolution, their mean value and their 
standard deviation, need to be known. 

This chapter will explain the nature of the wavefront distortions, starting with a 
short analysis of atmospheric turbulence. A more detailed analysis can be found 
in Roddier's Adaptive Optics in Astronomy [18] and Beckers' Adaptive Optics 
for Astronomy [T]. Methods for quantifying seeing are described in Chapter |4] 



Changes in the refractive index of air are essentially proportional to changes in 
the air temperature. Temperature inhomogeneities are produced when layers 
of different temperatures are mixed due to wind shears. The statistics of re- 
fractive index inhomogeneities follow the inhomogeneities of temperature which 
are described by Kolmogorov [12]. The model assumes that energy is inserted 
at low frequencies on large scales which are characterized by the outer scale of 
turbulence. There is also an inner scale set by molecular friction. Energy is 
transported from the outer to the inner scale in a cascading way and is finally 
converted into heat [T]. 

The value of the refractive index n is of little interest when studying wave- 
front perturbations caused by the variation of n. The Kolmogorov model states 
that the variance of the difference between the refractive index in two different 
locations, separated by the three-dimensional separation vector p, is given by 



where r is a three-dimensional positions vector and p = \p\. Dn is called the 
index structure function and the approximation is only valid when p is smaller 
than the outer scale of turbulence. The value for this outer scale is highly 



2.1 Atmospheric turbulence 



D^ip) = (|n(r) - n{r + p)\') ^ Cl{z)p'/\ 



(2.1) 
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debated and quoted values span from a few tens of centimeters up to kilometers 



The power law of Eq. (2.1) has experimentally been found to be quite accurate 



over distances less than 1 meter. That is, Eq. (2.1) is certainly valid for small 



telescopes but it is likely to be more inaccurate for large telescopes [H 



C^(2;) in Eq. (2.1) is the index structure coefhcient and varies over distances 



at a much larger scale than over the inhomogeneities. The integral over Cff{z) 
along the line of sight gives a measure of the total wavefront distortion at ground 
level. 

The temporal variance, how fast the index fluctuates with time, is as important 
as the spatial variance for adaptive optics systems. This temporal variance can 
be described in the same way as the spatial, by a temporal structure function 
describing the variance of the difference between n at time t and a later time 
t + T. Assuming that variations in n live longer than the time it takes for wind 
driven inhomogeneities to cross over the telescope, Taylor expansion gives 



n{r, t + t) — n{r — vt, t), 



(2.2) 



where v is the wind velocity. Equation ( 2.1 ) can with this assumption be rewrit- 
ten as a temporal structure function 



Dn(vT) ^ {\n{r,t)^n{r,t + T)\^) ^ 

= {\n{r, t) - n{r - vr, t)\') = CI{z)\vt\^/^ 



(2.3) 



and it can be concluded that the temporal structure function can be obtained 



by simply substituting \vt\ for p in the spatial structure function Eq. (2.1 1. 



Wavefront distortions 

The effect of the temporal and spatial variation of the index of refraction in 
a turbulent atmosphere on the wavefront perturbation are of great interest for 
adaptive optics. As described in Sec. |1.1| the distortion is caused by varying 
optical path lengths through the atmosphere, due to the inhomogeneity in n. 
The optical path length is given by 

n{z)dz, (2.4) 

where n is integrated over the line of sight. The air refractive index is fairly 
wavelength independent in the range from visible to near infrared and hence also 
the optical path length. However, the wavefront phase fluctuation is wavelength 
dependent and given by 

V^k j n{z)dz = kS, (2.5) 

where k is the wave number which varies as the inverse of the wavelength, 
k = 2tt/X. The distortions that are corrected for are however the difference in 
path length and they can therefore be compensated for at all wavelengths (since 
the path length is wavelength independent). 
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The difference between two points at tlie telescope entrance aperture is (again) of 
more interest than the absolute wave-front phase. The variance of the difference 
in the structure function between the phase in two different locations, separated 
by the two-dimensional separation vector ^ is given by 



(2.6) 



where x is the two-dimensional position vector and ^ — \^\ is the distance 
between the two points at the telescope entrance aperture. Using Eq. (2.5) in 



Eq. (2.6) makes it possible to express the phase structure function in terms of 



the index structure function integrated along the line of sight. Using Eq. (2.1) 
and integrating gives an equation for D^p in terms of C^(z) 



oo 

D^ii) = e J Cl{z) Az J ie + z2)i/3 _ 



(2.7) 



where the integral / is evaluated as 

2r(i)m), 



r(i) 



givmg 



D^{£,)^2mk' j Cl{z)Az^^'\ 



(2., 



(2.9) 



(z) is only dependent on height since the atmosphere generally is considered 



to be stratified in plane parallel layers. Equation ( 2.9 ) can therefore be rewritten 
as 



D^{i) ^ 2.91 



cos(7) 



Cl{h) Ahfl^. 



(2.10) 



where I ) is the air mass and 7 is the angular distance from zenith. Equa- 



tion (2.101 is usually rewritten as 



D^{i) = 6.88(C/ro)5/^ 
where the Fried parameter |8], 



ro 



0.423- 



Cij{h) Ah 



-3/5 



cos(7) 

is a useful quantity when characterizing seeing quality. 



(2.11) 



(2.12) 



Equation (2.111 describes the spatial distribution of the wave-front distortions. 



Using Eq. (2.3), a temporal structure function for the phase can be found by 

(2.13) 



substituting vt for ^ in Eq. (2.11 ) 



D^{vt) = 6.88(^;T/ro)^/^ 
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where v is the velocity with which the wavefront-phase propagates. This velocity 
is a weighted average of the layers' velocities. The effect of the finite response 
time of an AO system can be calculated with Eq. (2.13). If a correction to the 
wavefront is applied at t + r when the phase is measured at t, the mean square 
phase error is given by 



,{t) = 6.88(wT/ro)^/3^ 



(2.14) 



Equation (2.14) can be used to calculate a maximum allowed time delay before 



applying correction 



ro = 6.88-t^=0.314':". 

V V 



(2.15) 



This delay is called the Greenwood time delay and the corresponding frequency, 
fc is called the Greenwood frequency pD] . 

The isoplanatic angle is also a frequently used quantity to describe atmospheric 
turbulence. It defines an angular distance 6 from the source that is used to 
sense the wavefront to objects around this source for which the compensation 
still will be good. The image quality decreases as the angular distance increases. 
The mean square error in the wavefront for a layer at a distance h/cos{'j) is 
obtained by replacing ^ with 9h/ cos(7) in Eq. (2.11). For this to be valid not 



only for one layer, h should be replaced by the weighted average h of the layer 



altitudes (similar to what was done with v in Eq. (2.14)) giving 



,{9) = 6.88 



0h 



ro cos(7) 



5/3 



(2.16) 



To have an anisoplanicity rms error less than 1 radian the isoplanatic angle must 
be less than 



On = 6.. 



.3 To cos(7) 



0.314 



ro cos(7) 



(2.17) 



2.3 Wavefront representation 

It is often useful to describe the wavefront for a circular aperture in terms of the 
orthogonal set of Zernike polynomials, Zj{p,9). They can be expressed either 
in Cartesian {x,y) or Polar {r,9) coordinates. The polar form is most common 
and will be used here. 

Zevenj = V^T+IR^ {r)V2 cosimO) 1 

Zoddj = V^^K{r)V2sm{me) / ™^ (2.18) 



where 



(n-m)/2 t T\s( \\ 

R^(r] ^ V ^ ' — r"-2'' (2 19) 

"^^ ^ s![(n + m)/2-s]![(n-TO)/2-s]! ^' ' 



s=0 
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A wavefront phase-distortion over a circular aperture of radius R can be ex- 
pressed as a Zernike polynomial expansion, 



ip{Rp,9) = ^ajZj{p,e), 



(2.20) 



where p = r/R and aj some coefficients (for more details see [16]). 

The polynomials are set up so that the first few directly relates to optical phe- 
nomena such as tip, tilt, defocus, astigmatism etc. The first eight Zernike poly- 



nomials are described in Table 2.1 together with their optical error. Figure 2.1 
shows visualisations of the first eight modes. 



j 


Zernike mode Zj 


Aj (rad^) 


Optical aberration 


1 


1 


1.0299 


Piston 


2 
3 
4 


2r cos{9) 
2r sm{e) 
V3{2r^ - 1) 


0.582 
0.134 
0.111 


Tip/tilt (X-axis) 
Tip/tilt (Y-axis) 
Defocus 


5 


VEr^ sin(26') 


0.0880 


Astigmatism (±45°) 


6 


y6r2 cos(26l) 


0.0648 


Astigmatism (0790°) 


7 


78(37-2 _ sm{e) 
V8(3r2 _ cos{9) 


0.0587 


Coma (Y-axis) 


8 


0.0525 


Coma (X-axis) 



Table 2.1: The eight first Zernike polynomials. The constant Aj is used to determine 



the mean square residual error after correcting j Zernike modes, see Eq. (2.221. The 
last column describes the optical aberration that is associated with mode j. Table 
adapted from |16| . 



Figure 2.1: Visualisation of the first eight Zernike modes Zj. From left to right: 
Piston {Zi), Tip/tilt {Z2 and Z3), Defocus {Z4), Astigmatism (Z5 and Zq), Coma (Z7 



and Zs). Picture adapted from http://en.academic.ru/dic.nsf/enwiki/973636 



The correction for the first N modes can be written as 

N 



:{Rp,9) = Y,ajZ,{p,i 



(2.21) 



Fried [5] and Noll TF] found that if this correction is applied to a distorted 
wavefront, and assuming that Eq. (2.111 holds, the mean square residual error 
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can be defined as 

N /n\5/3 



-| = (^^)-E(KI')=^.\^)^ (2-22) 



where D is the aperture diameter and ro the Fried parameter from Eq. (2.121. 



The first eight values of Aj are found in Table 2.1 Noll [TB] found that for 



iV > 10 the the mean square residual error is given by 



cr| w 0.29447V" ( — j . (2.23) 



5/3 



The residual error without correction (af) is given by 
/ r)\ 5/3 

cr? = 1.0299 f—j , (2.24) 

so the rms phase distortions will be about 1 radian if a telescope aperture of 
size ro is used. 

Another method for decomposing wavefront errors is the Karhunen-Loeve ex- 
pansion. Like Zernike polynomials, Karhunen-Loeve (K-L) modes are orthogo- 
nal. In addition, they are statistically independent modes which are numerical 
instead of analytical (Zernike modes) and which can be approximated as sums 
of Zernike modes by diagonalization of the Zernike covariance matrix. 

The difference between a Zernike mode and its related K-L mode becomes larger 
with increasing order. The mean square residual error for an AO system com- 
pensating for the first TV K-L modes decreases faster than for a system compen- 
sating for Zernike modes. This means that it is more efficient to compensate 
K-L modes than Zernike modes. 
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rd Chapter 

Adaptive Optics 



Adaptive Optics (AO) is real time compensation of the image degradation, that 
is a technique to correct distorted wavefronts in real time. The technique was 
first proposed in 1953 by Babcock and later in 1957 by Linnick (independently). 
An adaptive optics system uses a deformable mirror to instantaneously correct 
the wavefront distortions [l8] . 

Figure [XT] schematically describes an adaptive optics system. 




Figure 3.1: A sketch describing a typical adaptive optics system using a deformable 
mirror as wavefront corrector and a Shack-Hartmann wavefront sensor to detect the 
distortions. The wavefront sensor sends signals to the deformable mirror forming a 
closed-feedback loop. Most light is transmitted by a beam splitter into the science 
camera while some is transmitted into the wavefront sensor. 



This chapter will describe the main parts of an adaptive optics a system but 
the first section will clear out the confusion between adaptive and active optics. 
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3.1 Active Optics 

The difference between adaptive and active optics can be a bit confusing. Active 
Optics is commonly used to describe ways of controlling wavefront distortions 
that are caused by optical, thermal or mechanical effects in the telescope itself. 
Active Optics is slow compared to Adaptive Optics. The effects that the active 
optics system is built to compensate for vary on a long time scale compared to 
seeing which changes rapidly. 

Active Optics works at bandwidths less than 1 Hz while Adaptive Optics uses 
bandwidths from 10 to 1000 Hz. Active Optics can use large primary mirrors 
for wavefront correction whereas Adaptive Optics need to use small mirrors 
optically conjugated to a lenslet array. [T] 



3.2 Wavefront sensor 

The wavefront sensor (WFS) is one of the basic elements in an adaptive optics 
system. The wavefront phase aberration to be corrected must first be mea- 
sured. This is done by a wavefront sensor. The wavefront has to be sensed with 
enough spatial resolution and enough speed for real time compensation of the 
atmospheric seeing. 

No wavefront-phase sensors for visible light exist today [TS]. The deformation 
is measured indirectly. This can be done in various ways. 



3.2.1 Shack-Hartmann wavefront sensor 

The Shack-Hartmann^ (SH) wavefront sensors are the most commonly used in 
astronomy [1 . A SH WFS consists of an array of small lenses (subapertures) 
in the pupil plane. The subapertures (SA) image only a part of the wavefront 
so that they only sense the local wavefront tilt. Each SA forms an image of 
the source at its focus, if the wavefront is plane. If the wavefront is disturbed, 
each SA receives a tilted wavefront and forms an off-axis image in its focal plane. 
The measurement of this displacement is a direct measure of the local wavefront 
tilt. By using an array of subapertures the wavefront slope can be obtained at 
a fixed number of points and the shape of the wavefront can be reconstructed. 



Figure 3.2 shows the principle of a Shack-Hartmann wavefront sensor. 

There are several methods for measuring the positions of the images formed 
by the lenslet array. Calculating the center of gravity for each subaperture, by 
tracking the peak of the image with a artificial (or natural) point source, is one 
method that works for point- like objects. [15] . 



^The design of the sensor was developed in 1900 by Johannes Franz Hartmann and was 
modified in late 1960s by Roland Shack and Ben C. Piatt who replaced the opaque screen with 
a lenslet array. The wavefront sensor is therefore also known as Hartmann-Shack wavefront 
sensor. 
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Figure 3.2: Principle of a Shack-Hartmann wavefront sensor. The perturbed wave- 
fronts arrives at top. The picture shows a 1-dimensional cut through the subaperture 
array with 4 individual subapertures (SA). Each of these SA samples a different part 
of the incoming wavefront, sensing only the local slope. The slopes are drawn with 
dashed lines. The image displacement is shown in the leftmost subaperture. The 
image would have been centered if the wavefront had been plane, it is now offset by 
Ax. 



For extended objects, like the Sun, one has to measure relative image motion by 
comparing the same features in different subpupil images. Subaperture images 
are cross-correlated with a chosen reference image for various displacements. 
The displacement that yields the highest cross-correlation is used as the sub- 
aperture offset. 

Lofdahl [M] has evaluated the inherent accuracy of the most commonly used 
methods. 



3.2.2 Wide Field Wavefront Sensor 

The Wide Field Wavefront Sensor (WFWFS) is a SH with a large field of view. 
Because of this large FOV the seeing distribution from different layers (see 



Fig. 1.2) in the atmosphere can be tomographically mapped using the angular 
dependence between the different wavefront slopes. A wide FOV lowers the 
height for the first seeing layer but a narrow FOV gives better height resolution. 



The height variation is explained further in Sec. 4.2.2 



3.3 Wavefront correctors 

When the wavefront distortion is known it is up to the optical components of 
the AO system to correct it. The optical components are usually deformable 
mirrors. They are designed so that their shape can be adjusted in order to match 
the incoming wavefront distortion. The function of the optical components of 
an AO system is often divided among two components, a tip-tilt mirror and an 
adaptive mirror that corrects the higher order wavefront distortions 

The tip-tilt mirror is usually the first corrector. It is a flat mirror that can 
rotate on two perpendicular axes. The tip-tilt mirror only corrects the tip and 
tilt modes while the deformable mirror compensates the high-order aberrations. 
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Most of the atmospheric wavefront aberrations are located in the tip and tilt 



modes (see Table 2.1 ) and the tip-tilt mirror will therefore need a larger stroke 



than the deformable mirror. 

The deformable mirror is a high-order corrector and typically has less stroke 
than the tip-tilt mirror. The efficiency of the AO system will be reduced if this 
mirror has to correct for low order modes as well as high order modes since 
there will be less stroke left for the high-order modes. 



Multi-conjugate adaptive optics 

As described above, a conventional adaptive optics system consists of a de- 
formable mirror and a tip-tilt mirror integrated in the telescope's optical path, 
a wavefront sensor and a control computer. 

This can be extended to a Multi-conjugate Adaptive Optics system by inserting 
several deformable mirrors as well as wavefront sensors. The deformable mirrors 
in the MCAO system are optically conjugated at different altitudes. 

The Sun is well suited for MCAO since it is an extended object and it displays 
small scale structure everywhere on its surface. A "natural guide star" can be 
found wherever the telescope points. ^Sj 
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th Chapter 

Methods for 
quantifying seeing 



Several methods for quantifying seeing exist. Differential image motion monitors 
and shadow band ranging will be described in this chapter. 

The chapter ends with a description of the optical setup at the SST. 



Differential image motion monitor 

The differential image motion monitor is the most widely used method for quan- 
tifying seeing. Wavefront slope differences are measured over two small pupils 
some distance apart. The differential method makes it insensitive to tracking 
errors. It has been shown that the differential motion exceeds the absolute 
motion when the distance between the two apertures equals a few times their 
diameter [THj. 

Sarazin and Roddier |19j presented the following approximate equations for the 
variance of differential image displacements measured with two subapertures 
with a diameter D and separated, along the x-axis, by a distance s (nomencla- 



ture is chosen to match Sec. 4.2 and Fig. 4.1 ) 



af = {{x{s) - x{0)f) = 0.358X%'/'d-^/^ (^1 - 0.541 (^) , (4.1) 

for longitudinal displacement (a long a subaperture row) and where tq is Fried's 
parameter. The transverse displacements (perpendicular to the subaperture 
row, a long a column) are given by 



= {{y{s) - y(O)r) = 0.358X\^'/' D-'^' ( 1 - 0.811 (-^) 



-1/3 



(4.2) 



Using the notation of Fried [S| Equation (4.1 ) and Eq. (4.2) can be written 



{{x{s) ^ x{0)f) = 0.358X''r~^^'^D-^/^I{s/D,0), 
{{y{s) ~ y{0) f) = 0.358AV^'^"'/'^(^/A|), 



(4.3) 
(4.4) 
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where I{s/D,0) and I{s/D, |) are tabulated in Fried's paper [5]. The function 
/ is normalized so that it approaches unity when s approaches infinity. It is also 
symmetric so that Iis/D,0) = I{-s/D,0) and I{s/D, f ) = I{-s/D, f ) 



4.2 S-DIMM AND S-DIMM+ 

The Solar DIMM (S-DIMM) was first used in the Yunnan Observatory and uses 
differential motions of the solar limb measured with two small subapertures 
[15] ■ The S-DIMM+j which can be seen as a natural extension of the DIMM 
and S-DIMM, includes measurements that are sensitive to height variation of 
seeing. The calculations made by Fried must be extended to account for the 
averaging effect of using a large field of view for wavefront sensing and to get 
more accurate estimates of / when s/D approaches zero |20| . 



4.2.1 Extension of Fried's calculations 

Scharmer and van Werkhoven j2Qj extended the calculations of Fried [9] to 
include s/D < 1 and s/D as large as 50, the latter corresponding to a distance of 
5 m for subapertures with a diameter of approximately 10 cm. Their calculations 
agree perfectly with those of Fried for s/_D = 1 but they deviate with increasing 
s/D and differ by 10 % for s/D = 10. This suggest uncertainties in Fried's 
numerical integration and was suspected already by Sarazin and Roddier |19j . 
For s/D > 7.5, the calculations of Fried show a small decrease of longitudinal 
differential image displacements. Scharmer and van Werkhovens's calculations 
behaved similar when the step size of one of the integrating variables were to 
large and they therefore decreased that step size with approximately a factor 
of 80 compared to the calculations that essentially reproduced the results of 
Fried. After that improvement the calculations showed monotonous increase of 
I{s/D,0) and I{s/D, ^) with s as long as s/D < 300. That is a much larger 
range than needed for their analysis and it could also be concluded that their 
calculations must be more accurate than Fried's. 

Scharmer and van Werkhoven conclude that the theory developed by Fried can 
be used for the S-DIMM-I- by replacing the separation s with s+ah and the pupil 
diameter D with an effective diameter D^g- An approximately round binary 
mask is applied when measuring image displacements from the granulation with 
cross-correlation techniques in order to make the averaging area roundish. 



Figure 4.1 describes the relation between the wavefront sensor geometry and 
the two differential wavefront slope measurements at a certain height h. It is 
assumed that the first subaperture is located in the origin with a separation 
of s to the second subaperture. The field angle of the first measurement is 
furthermore assumed to be zero while that of the second measurement is a 
and that the tilt only is along the x-component of the image displacement. The 
measured quantity 5xi for the first measurement, corresponding to contributions 
added from the A^ layers located at heights is 

N 

(5x1(5, 0) = ^(x„(s)-x„(0)), (4.5) 
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while the quantity for the second measurement will be 

N 

6x2(3, a) = ^(a;,i(s + ahn) - a;„(a/i„)). 

71=1 

The covariance between dxi and 6x2 is given by 

N 

{6x16x2) = ^((a;„(s) - Xn{0)){Xnis + ahn) ~ Xniahn))), 



(4.6) 



(4.7) 



because of the assumed independence of contributions from different layers. The 
average brackets, (...), in Eq. (4.7) denote averages over many exposed frames. 
By evaluating the four terms in Eq. (4.7) one by one and by taking advantage 
of the assumed homogeneity of the statistical averages at each height /i„, this 
equation can be rewritten as a combination of three variances of differential 
image displacements, 

N 



{6x16x2) = '^{{xn{ah„ - s) - a;„(0))^)/2 - 



{{xniahn + s)- x„(0))2)/2 - ((x„(a/i„) - Xn{Ojf). (4.8) 

When a = 0, the last term disappears and the first two terms are equal which 
means that the equation describes a conventional DIMM, except that the contri- 
butions from the high layers are reduced by the averaging of image displacements 
from the relatively large FOV used with the S-DIMM-I- method. 



By combining Eqs. (4.3 4.4| and (4.8), the following equations can be obtained: 

N 



{6x16x2) = ^ CnF^{s, a, /i„), 

N 

{6yi6y2) = ^ c„i^j,(s, a, /i„). 



(4.9) 
(4.10) 



where 



F,{s,a,hn) - /((a/i„ - s)/-Deff, 0)/2 + 

I{{ahn + s)/D,s, 0)/2 - I{{ahn)/D,s, 0), (4.11) 
Fy{s,a,hn) = / ((a/i„ - s)/L'cff, I) /2 + 

/ (^{aK + s)/D,ff, I) /2 - (^{aK)/D,ff, |+) , (4.12) 



and the coefficients c„ come from the original DIMM equations Eqs. (4.3 4.4), 
rewritten to match the S-DIMM+ method. 



0.358 



A 



5/3 



A 

F>cs{hn) 



1/3 



(4.13) 
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4.2.2 Height variation 



Figure 4.1 describes the relation between the WFWFS geometry and the contri- 
bution to the image displacement measurements from a specific height. Large 
separation between the subapertures gives better height resolution but smaller 
separation gives a higher maximum height. A large FOV lowers the minimum 
height but a small FOV, on the other hand, gives better height resolution. 

The size of the subapertures is also important. Smaller subapertures give better 
height resolution, but it is limited by the FOV for large heights and bigger 
subapertures lowers the noise. 




Figure 4.1: The relation between the wavefront sensor geometry and the contribu- 
tions to the two differential image displacement measurements from a height h. s is 
the separation between the two subapertures and a is the measured relative field angle 
between the two subfields. The beams diverge with height due to the non-zero FOV. 
Picture adapted from |20) . 



Eqs. (4.9 4.10) can also be expressed in terms of the turbulent strength of 
each layer, dh. The atmospheric structure constant Cff can be found by 



rearranging Eq. (2.12) into 



/ X \ 5/3 

C^d/i = 0.0599 cos(7) f — j X^^^, 



(4.14) 



with which Eq. (4.13) can be rewritten as 



5.98 / 1 



cos(7) \Dcff{K 



1/3 



Cj^dh. 



(4.15) 



Scharmer and van Werkhoven found that the S-DIMM-I- (for the SST) should be 
able to distinguish seeing at the pupil where h — and at a height of /i = 500 m 
above the pupil and that the angular resolution is enough to allow measurements 
up to a height of approximately 20-30 km. 
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Shadow band ranging 

As mentioned in Sec. |l.l| stars seem to twinkle when observed with your naked 
eye. The proper name for this twinkhng is astronomical scintillation. The light 
from a star is focused as a point at the back of your eye. As shown in Fig. |4.1| 
the beams gradually widens closer to the star (due to a non zero FOV). The 
size of this beam is still very small at the height of the turbulent layers and the 
star seems to twinkle since the light is disturbed. 

Seykora [53] showed that the scintillation of solar light was closely correlated 
to seeing and measuring this scintillation using a photodiode is hence another 
method for quantifying seeing PJ [23] . This is done by a scintillometer which 
detects and evaluates the intensity fluctuations of the transmitted signal. 

By assembling several detectors in a scintillometer array, information about the 
structure function can be obtained |4]. This was further developed into 
an instrument called SHadow BAnd Ranger (SHABAR) which consists of a 
non-redundant linear array of scintillometers . 

A SHABAR measures the scintillation index, erf, which is defined as the variance 
of the relative irradiance fluctuations. The covariance between the SHABARs 
scintillometer units are then inverted to achieve the profile. 

Higher layer contribute less to he scintillation index and the SHABAR has 
therefore a limited range. A longer baseline for the scintillometer array is needed 
to achieve sensitivity on higher heights. 

The short baseline SHABAR that is mounted on the SST has its 6 scintillometer 



cells inserted in a 50 cm long bar shown in Fig. 4.2 |24] 
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4.4 Optical setup at the SST 



Figure 4.3 describes the present optical setup used at the SST. 



Field 

WFWFS stop 




Wavcfront Sensor 



Blue 

camera 

tower 



Beamsplitter 



Optical setup at the SST 



Figure 4.3: The present optical setup for the SST. The incoming light is folded from 
almost vertical to horizontal when passing the tip/tilt (TT) mirror and the deformable 
mirror (DM). The light is split at ~ 500 nm into two different science beams. A 
correlation tracker (CT) is located in the blue beam while CRISP and the WFS are 
located in the red. The WFWFS is mounted directly under the vacuum system so 
that the WFWFS beam doesn't pass the AO-system. Drawing adapted from |26| . 



The telescope pupil is re-imaged onto a deformable mirror (DM) via a field lens 
(LO) and a tip-tilt mirror (TT). The tip-tilt mirror deflects the beam by 60° 
and the deformable mirror deflects it by 30°, that is the beam is folded from 
nearly vertical to horizontal. The incident angle of 15° on the adaptive mirror 
makes the pupil image nearly circular (cos(15°)=0.97). A re-imaging lens (LI) 
magnifies the image ~ 2.2 times, making an image scale of ~ 0.041 arc sec per 
pixel on the CCD in the blue beam and ~ 0.059 arc sec per pixel in the red 
beam. 

The electrode pattern of the 37-electrode bimorph deformable mirror is placed 
upon each of two thin plates of piezo electric material which are glued to the 
mirror. When voltage is applied to an electrode, one of the plates will expand 
while another will contract, bending the bimorph. |26j 

The horizontal part of the beam is split by a dichroic beam splitter into two 
different science beams. A correlation tracker is located in the blue beam and 
the SH WFS is located in the red beam (see Fig. 4.3). The correlation tracker 



controls the tip-tilt mirror and the WFS controls the differential mirror. 

Another beam splitter (partially polarizing) splits the light so that 8 % goes 
into the WFS. The wavelength range is restricted by a 100 nm wide interference 
filter that is centered at 550 nm. This filter will also improve the image quality 
by preventing chromatic errors. The field of view is defined by an adjustable 
field stop at the secondary focus (FS2) so that the subaperture images don't 



20 



4.4. OPTICAL SETUP AT THE SST 



overlap. The collimator lens (L2) that is located after the secondary focus, but 
before the subaperture array, produces a parallel beam and re-images the DM 
onto the subaperture array. The subaperture array re-image the secondary focus 
onto the detector. [26] 



4.4.1 WFWFS DESCRIPTION 

The WFWFS is mounted immediately below the vacuum system. The field 
of view of the WFWFS is placed outside the science field of view and light is 
deflected horizontally so that the WFWFS beam does not pass the SST AO- 
system. The optics of the wavefront sensor consists of a field stop, a collimator 
lens and an array with 85 hexagonal subapertures within the 98 cm diameter 
pupil. This layout is shown in Fig. |4.4a| The diameter of the subapertures is 
equivalent to 9.8 cm. The WFWFS has a maximum field angle of 46.4 arc sec. 
The image scale is 0.344 arc sec per pixel. 




(a) 
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Figure 4.4: (a) Layout of the 85 fully illuminated hexagonal subapertures. The circle 
represents the 98 cm SST aperture. The two subapertures Sri and Sr2 corresponds to 
selected reference subimages with high RMS contrast. Arrows point to the subaper- 
tures s — and s for which differential image shifts are measured, (b) Corresponding 
subimages with masks (not to scale), indicate the two subfields at field angles separated 
by Q. Figures adapted from pU] . 
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Figure 4.4a shows the 85 hexagonal subapertures (SA) and Fig. 4.4b show the 
corresponding subimages. Subimages are located on the CCD while subaper- 
tures are located in the aperture plane (on the lenslet array). Every subimage 
is divided into several subfields. This is shown in Fig. |4.4| 
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Data processing 



After observing with the telescope, the images taken with the WFWFS are in- 
spected to see if they are good so that they can be used for further analysis. If 
not, a warning or even an error is triggered. This inspection is done automat- 
ically directly after a data set is taken but can of course be started manually 
afterwards. The different data sets are then reduced using a program called 
Astooki which, after several steps, returns a S-DIMM-I- covariance map. 

This chapter describes the inspection and some statistics made out of the inspec- 
tion as well as the different reduction steps. The chapter includes a description 
of how data from the SHABAR (also mounted on the SST) were processed be- 
fore they were compared with the WFWFS data and ends with a description of 
changes in the reduction. 



Autoproxy is an automated proxy software, written by T.I.M. van Werkhoven, 
which acquires, checks and stores WFWFS data [28] . The checking or inspection 
starts after the data has been acquired. The quality of the data is checked and 
some statistics on the frames are calculated, for example min, max and average 
intensity, RMS, saturation etc. Some global statistics, such as the standard 
deviation of the frame statistics, are also computed. 

The values of theses parameters are used to judge the quality of each frame and 
assign a status code. Warnings are indicated by positive status codes and errors 
give negative codes. Status code implies a 1700^ frame. The threshold values 
for each status code can be found in Table 15.11 

A log file is written in the end of the inspection including the statistics of each 
frame as well as the assigned status code. The log file also states the first, last 
and best frame, together with information about the number of frames that 
where classified as good, warning or bad. This can then be used to classify the 
whole data set. 

Frames with errors, which during the inspection were marked as bad, are re- 
moved by Autoproxy to save storage space. This setting can of course be 
changed when, for some reason, also bad data need to be stored. 



Inspection 
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G4- n-i-nc- 

Status 
code 


Criteria 


Image 
Warn Bad 


Flats 
Warn Bad 


Darks 
Warn Bad 


± 4 


Too many saturated pixels 


0.05 % 


0.2 % 


0.1 % 


0.2 % 


0.01 % 0.1 % 




Too low max. intensity 


<2000 


<800 


<2000 


<800 




± 16 


Intensity problems, average 
of frame 


>3200 
<2000 


>3600 
<600 


>2800 
<1000 


>3200 
<400 


>150 >500 


^± 32 


Intensity problems, average 


>3200 


>3600 


>3200 


>3600 






of SA 


<2000 


<600 


<1600 


<600 




± 64 


Frame avg. > SA avg. 


70 % 


64 % 


70 % 


64 % 




^ +128 


Frame after series of bad 
frames, probably also bad I 












± 256 


Too large SA avg. deviation 


>80 


>200 








^± 512 


Too high RMS 


6 % 


10 % 


2 % 


4% 





Table 5.1: Criteria for file inspection. The plus and minus sign indicate that a frame 
has a warning status (+) or a bad status (— ). If a frame (image, flat or dark) has 
saturated pixels it will be assigned the status 4 (warning or bad depending on the 
amount of saturated pixels). Status 8, 16 and 32 are applied to frames if they have 
too low or too high intensity (maximum or average) . Status 64 is assigned if the frame 
has an average intensity that is higher than 64 % (warning) or 70 % (bad) of the 
average intensity of the subaperture. Status code 128 is only assigned as a warning 
to frames that follow long series of bad frames. Status 256 is applied if the difference 
between the frames subaperture (SA) value and the average value for all SA in that 
set is too big. Status 512 is assigned if the relative RMS value of the pixels within the 
subapertures are too large. A dash (-) indicates that the criterium is not available for 
that file type. 



During this master thesis project, several scripts were written in order to make 
it easier to inspect the result of the automated inspection. A script that checks 
the number of bad frames for every status code is shown in Appcndix |A.l[ There 
is a corresponding script for warnings. 

5.1.1 Status codes 

The 8 available status codes are described in Table [5T] with their corresponding 
thresholds. A frame can be assigned with one of these status codes or with a 
combination of them. Status codes are " combined" with the binary OR operator 
in order to get unique status codes. Status code 48 for example is assigned to 
a file that has a waring due to to low or to high average intensity in the frame 
(status 16) and in the subaperture (status 32), that is 16|32 = 48. 

The combination of status codes is unique as long as the frame only has warn- 
ings. When the frames start to get errors, some status codes can be a "combina- 
tion" of different statuses. Status —200 for example can either be 16|32| — 256|8 
or 16|32|256| — 512|8. A frame with one error status code will though always 
have a resulting code that is negative (since e.g. 512| — 8 = —8). 
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5.1.2 Change in inspection 

A lot of frames from 2009 were browsed through with the ANA browser and 
SST Java Browser available on the local computers to see if the thresholds 
for the inspection seemed to have the right values. During this browsing, it was 
found that some of the frames marked as good were shifted or had overlapping 
subapertures. 

Inserting more field stops in the optical setup and making sure that the WFWFS 
does not take data when something is blocking its field of view would decrease 
the number of such frames. However, the remaining bad frames still need to be 
found without manually inspect a million files. 

A small change in the threshold for one of the already existing status codes 
(64) marks these frames as warning or bad frames depending on the amount of 
badness. 

Status code —64 (indicating a bad frame) originally had following threshold 
frame average > subaperture average. This was changed to the more strict 
thresholds frame average > 0.7 * subaperture average resulting in status 
code —64 (bad) and frame average > 0.64 * subaperture average result- 
ing in +64 (warning). 



5.1.3 Intensity and RMS 

Frames are assigned with a warning (or an error) if they have too low intensity. 
Since quite many frames were assigned with this type of warning, a script was 
made to check which intensities these frames actually had, that is if they were 
close to the given thresholds or not. Some frames were then examined by eye. 

Frames which have a too large relative RMS intensity for the pixels within a 
subaperture are assigned with a warning (or an error). A similar script was 
written to check if some frames where deleted because of high RMS and if they 
where close to the threshold. 



5.2 Reduction 

The reduction of WFWFS data is done in several steps using Astooki - the 
Astronomical Toolkit^ which is a Python program made to process astronomical 
data [37|. The steps are explained below. 

During this master thesis project a reduction script was written in order to make 
all of the following steps in one go by calling the different tools of Astooki. The 
reduction script loops over all data sets for all specified days so that data from 
a whole year can be reduced in one go if wanted to. 

At first the script tries to find the fiat field and dark frame that are closest 
in time to a particular data set. The fiats and darks, together with the image 
frames, are needed as input for some of the Astooki tools. 

Code snippets from the reduction script can be found in Appendix [B| 



25 



CHAPTER 5. DATA PROCESSING 



5.2.1 Initial setup 

Making a subimage mask 

The first step is to create a subimage mask. A set of coordinates for where the 
subimage is located are calculated. The units are arbitrary but they should be 
consistent, pixels (on the CCD) were used for this mask. 

The radius (1024 px, with a 2k by 2k pixel CCD) of the pattern is specified 
as well as the pitch (194 px, 166 px) of subimage positions. The centroid co- 
ordinate, that is the coordinate in the middle of the subimage, and the origin 
coordinate (the coordinate in the lower left corner of the subimage) are cal- 
culated for every subimage. In order to match the subimage pattern with the 
corresponding hexagonal subaperture pattern, an offset is specified so that every 
odd row is shifted a half subimage in x-direction. 

It is checked if the obtained coordinate (with offset) lies within the aperture 
bounds. The telescope aperture is circular, and since the subimages on the 
CCD correspond to subapertures in the telescope aperture, the subimages also 
have to fit within the circular bound. 

The coordinates (with applied offset) for the subimages that where located 
within the circular bound are saved in two arrays of subimages (one with origin 
positions and one with centroid positions) that describes the subimage pattern. 



Part of the resulting pattern is shown is Fig. 4.4b 

The subimage mask is made with equal parameters for all data sets and days. 
The mask will be optimized for a particular data set in next step. 



Optimizing the subimage mask 

After the "standard" subimage mask is generated, it needs to be optimized since 
it is possible that even a mask created with the right parameters won't match 
the subimages perfectly. In this step, the mask is matched to a fiat field image 
(the one closest to that data set). 

Figure |5.1| explains how this fitting is done. The range, in which the vertical 
and horizontal slices are taken, is 1/2 subfield larger, in every direction, than 
the compared subfield. This is to make sure that the range covers the whole flat 
fielded image. Two slices are taken out within this range, one in the x-direction 



and one in the y-direction (see Fig. 5.1a I 



In order to quickly find the borders between the subimages, these slices are 
averaged (the x-slice over the height and the y-slice over the width) to one pixel 
as shown in Fig. |5.1b[ The routine defines the border as the first pixel to left 
and right from the center of the x-slice with an intensity less than 80% of the 
maximum intensity. This procedure gives the width of the subimage. The height 
of the subimage is then given by the first pixel to the top and bottom from the 
center of the y-slice that is less than 80% of the maximum intensity. 

Optimizing is done for all subimages. All sizes are stored and the mean subimage 
size is calculated as well as the standard deviation. 

This routine is slightly sensitive to specks of dust on the fiat field. 
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(a) (b) 

Figure 5.1: (a) Slices in the optimized subimage mask. The red square indicates the 
range in which the two slices are taken out. The green rectangular area shows the 
horizontal slice and the blue rectangular area shows the vertical slice, (b) Slices are 
averaged down to one pixel to faster find the borders. 

5.2.2 Calibration 

The reference position of each subimage needs to be known in order to compare 
the corresponding subfields in different subimages. It can not be assumed that 
pixel (x,y) in subimage A corresponds to the same position in subimage B (due 
to static aberrations). The problem can also be addressed the other way around: 
given a granule on the Sun, it has to be known at what pixel that granule is 
located in each of the subimages. 

The "static offsets" are measured by taking a large field of view (almost the 
complete subimage) in one reference subimage and then cross-correlate that 
with all subimages. That results in N shift vectors, where N is the number of 
subapertures. It is possible to use multiple subapertures as references in order to 
get better results. This should provide the same data and will give an indication 
of the reability of the shift measurement. [57| 

Making a big subfield mask 

The subfield mask is a grid of coordinates relative to the subimage and subfields 
cannot lie outside the subimage. The size of the subimage, which was calculated 
in the step before, is therefore used as an input parameter. 

A border of 30 x 30 pixels is also supplied in order to keep a guard range free at 
the edge of the subimage. This is needed in next step when the cross-correlation 
is done. 

The effective size of the big subfield mask will be the size of the area inside the 
border (30 x 30 px) that is applied to the optimized subimage size (from the 
step before). 
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Measure static offset shifts 



This part of Astooki, that calculates the image shifts between different subfields 
and subimages, is implemented in C. 

The big subimage mask and the subfield mask are, together with flat flelded and 
dark corrected data, used as input. The four subimages with the highest RMS 
value are selected as reference subimages. Each reference is compared with all 
subimages (even with itself) and a correlation map is calculated for every subim- 
age. The cross-correlation is done with the absolute difference squared method, 
Eq. (5.1 1, and the subpixel maximum of the correlation map is calculated with 
a 9 point interpolation method. 



Correlation map = - X! 1^^^' ~ 3ref{x + i,y + j)\\ , (5.1) 



where i and j is the shift (in x and y direction) and ^rcf is one of the reference 
images. The minus sign is to ensure a maximum for best match. 

The image shift is calculated for all subimages in all frames and stored in a 
quite big 5 dimensional matrix of data (^frames, preferences, #subimages, 
#subfields, #directions). 

The differential image shifts that are calculated with cross-correlation techniques 
suffer from a lack of absolute zero-point reference. It is in practice impossible 
to define such an absolute reference precisely, which is furthermore not needed 
since one can rely on the differential image shifts to average to zero over a 
sufficiently long time interval. |2D] 

One data set consists of 1000 frames, which corresponds to approximately 110 
seconds of WFWFS data. It is assumed that the seeing induced differential 
image shift averaged over these 110 seconds is zero. The average shift measured 
is therefore subtracted from the 1000 frames of the burst. 

The variance and the mean of these 1000 frames are calculated. The mean is 
the shift offset used in the next step (x and y offset for all subimages, averaged 
over frames). 



Updating subimage mask with calculated offset 

The previously calculated static offsets are now used to correct the subimage 
mask (that was calculated earlier). The origin position for every subimage 
mask is subtracted by the offset (that was obtained in the previous step) for 
that subimage. 

The maximum offset for all subimages (in one frame) is calculated and added to 
every subimage mask position. All subimage mask sizes are cropped with twice 
the maximum offset and the value for the smallest subimage size is stored. 

Each subimage mask will now point at the same location on the Sun. 
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5.2.3 Measure actual subfield/subimage shifts 

The subimage mask is now corrected with the static offsets and can be sub- 
divided in different subflelds so that the actual shifts for the data set can be 
calculated. 



Generate a subfield mask with smaller subfields 



This mask is created with the same routine as the previous big subfield but 
the resulting subfield mask will consist of several small subfields (instead of one 
big). 

Subfields can not lie outside the subimage and the smallest subimage size (found 
in the previous step) is used as input so that the subimage mask can be used 
for all frames. A 7 pixel wide border is applied to the subimage to keep a guard 
range at the edges. The effective size of the subimage mask (inside the border) 
is divided into subfields with a size of 16 x 16 pixels. The subfields overlap with 
50 %. Fig. 5.2 shows how the subfields are related to the subimage. 



Subimage size 

















Effective 
SI size 


SFl SF3 
SF2 



Figure 5.2: Upper left corner of a subimage mask with three subfields indicated. 
The border is 7 pixels wide and the subfields are located within the effective subimage 
size. The subfields are 16 x 16 pixels and overlap with 50 %. 



Measure subfield/subimage siiifts 

The optimized and updated subimage mask and the 16 x 16-subfield mask are 
used to calculate the image shift for each subfield in each subimage. The calcu- 
lations are done in the same way as for the static offset shifts described above, 
but with a circular mask applied to every subfield. 

Four subimages, the ones with highest RMS, are chosen as references. Every 
reference subimage, with all subfields, is compared with all subimages (and their 
subfields) creating a cross-correlation map with image shifts for all subfield- 
subimage pairs. 

The cross-correlation is done with the absolute difference squared method. 



Eq. (5.2) and with the circular mask applied. The subpixel maximum of the 



correlation map is calculated with a 9 point interpolation method. 

( y 

Correlation map = — j mask • \g{x, y) — g^dix + 2/ + j)| ] , (5-2) 
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where i and j is the shift (in x and y direction) and g^ef is one of the reference 
images. The minus sign is to ensure a maximum for best match. 

The image shift is calculated for all subimages in all frames and stored in a 
quite big 5 dimensional matrix of data (#frames, preferences, #subimages, 
^subfields, ^directions). 



5.2.4 S-DIMM+ ANALYSIS 

By using the calculated image shifts, the data can be fitted to a model in order 
to produce information about the actual seeing. This is done by the statistical 
S-DIMM+ method described in Sec. [4?2l 



Making a lenslet mask 

A subaperture mask of the lenslet coordinates is made with the same routine as 
the first generated subimage mask. Input parameters this time are the diameter 
of the microlenses and a radius that is slightly larger than the real aperture 
radius. Every subaperture has to fit 100 % within this radius and the radius 
chosen is therefore a little bit bigger. Some apertures in the real optical setup 
are a bit cropped. 

The centroid coordinate and the origin coordinate are calculated for every sub- 
aperture. In order to achieve a hexagonal subaperture pattern, an offset is 
specified so that every odd row is shifted a half subaperture in the x-direction. 
It is checked if the obtained coordinates lie within the circular aperture bound. 

The coordinates (with applied offset) for the subapertures that where located 
within the circular bound are saved in two arrays of subaperture positions (one 
with origin positions and one with centroid positions) that describe the sub- 
aperture pattern. This pattern is shown in Fig. |4.4a| 

The subaperture mask is made with equal parameters for all data sets and days. 
S-DIMM+ analysis 

This tool computes a covariance map for all subfield-subimage/subaperture 
pairs. Part of this Astooki tool is implemented in C. The previously calcu- 
lated image shifts are, together with the lenslet and subfield masks, used as 
input. 

The tool loops over all rows of subapertures (having the same y position). 
Within each subaperture, it loops over all subfield rows (subfields with the 
same y coordinate) and chooses a reference subfield. All other subapertures in 
that row, and their subfields are then compared to to the reference subaperture 
chosen. The comparison between the the subaperture-subfield pair found are 
compared as described in |4.2| 

The same procedure is then repeated for all columns. After this a S-DIMM-t- 
covariance map is produced. 
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5.2.5 Inversion in ANA 

The actual inversion is done in ANA and written by Goran Scharmer [2QJ . This 



step is not a part of Astooki. The unknown coefhcients c„ in Eqs. (4.9 4.10) are 
obtained by solving a linear least-square fit problem by minimizing the badness 
parameter given by 



N 



s,a \ n—1 

N \ 

+ [Cy{s,a)-J2<^nFy{s,a,h„)]^ W'(s,a), (5.3) 



where Cx{s,a) is the covariance matrix in absence of noise and the weight 
W{s, a) equals the number of independent measurements for each (s, a) |20| . 



As can be seen in Fig. 4.4a a large number of independent samples (subaper- 
tures) exist for a small separation s but only a small number exist for a large 
separation. This is of course similar for the small and large field angles a. To 
balance this varying number of measurements for different separation and an- 



gle, a weight {W{s,a)) is inserted in Eq. (5.3). This weighting is important 
especially for high-altitude seeing when the signatures are at small a and W is 
large [20] ■ 



Minimization of Eq. (5.3) leads to a linear matrix equation Ac = b for the 
coefficients c„. But this allows solutions with negative values of c„. A variable 
substitution 



(5.4) 



restricts the solutions to only positive values of c„. The problem is now a 
non-linear least-square problem, which is solved for the parameter (20j 

The height grids used for the inversion are found by calculating the inverse of 
matrix A, corresponding to a linear solution for c„. The grid is then chosen so 
that the noise sensitivity is minimized. 

The highest height possible with a FOV of 5.5 x 5.5 arc sec is 30 km and the node 
below that must be placed at 16 km in order to avoid high noise amplification. 
Better height resolution at large distances can only be improved with a smaller 
FOV. 

The heights used in this comparison are from 7 different layers (6 above the 



pupil plane). The correlation coefficient (Eq. (5.5)) was calculated to see how 
much these layers correlate with each other. The height grid used was h = 
[0, 500, 1300, 3200, 5800, 10800, 20000] m. 
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SHABAR DATA 

Data from the short-baseline SHABAR that is instaUed on the SST has been 
processed in order to make some comparisons with the WFWFS. 

The profiles measured by the SHABAR were integrated in order to obtain 
tq for a number of layers above h = 0. The boundaries between layers first used 
were hshabar = [1,30,60,90,120,200,500,1000,3000,8000,20000], giving ro for 
10 different SHABAR layers after integration. 

The correlation between these different layers was calculated using the linear 
correlation coefficient (r) also known as Pearson's r |17j . The correlation coef- 
ficient is given by 

J2{xi - x){yi - y) 

. (5.5) 



The value of r lies between -1 (complete negative correlation) and 1 (complete 
positive correlation), and indicates that there is no correlation at all between 
the variables (in this case no correlation between layers). 

The correlation matrix produced with these correlation coefficients shows the 
correlation between the different layers. Little correlation between the layers 
will result in a diagonal matrix. The location and number of boundary layers 
used for the integration of Cff were changed in order to get as little correlation 
as possible between the layers. 

In order to compare rg calculated from SHABAR measurements with tq calcu- 
lated from WFWFS following quantity was minimized 

roiWFWFS) -{jl^n- ro,niSHABAR)-^/^^ , (5.6) 

where rQ{WFWFS) is the Fried parameter calculated from WFWFS measure- 
ment at different heights and rQ^n{SHABAR) the Fried parameter for different 
SHABAR layers. The coefficients c„ are the unknowns and represent the best 
fit contributions from at different SHABAR layers used to fit WFWFS data 
at the pre-selected heights. 



Change of subfield size 

The size of the circular mask used together with the subfield mask when calcu- 
lating the image shifts was changed in order to improve the sensitivity at large 
heights. Better height resolution can be achieved with a smaller subfield size, 
but a larger subfield size on the other hand, gives lower noise and fewer failures. 

A circular mask with a diameter of 16 px, the same size as the subfield, was 
used in the first reduction. A second reduction was made with a 12 px mask and 
a third with a 8 px mask. The subfield size (16 x 16 pixels) was kept constant 
for easier comparison between the changes. Figure [573| shows how the different 
masks are related to the subfield size. 
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16 px 



12 px 



px 



Figure 5.3: The different "circular" masks that where applied to the subfields when 
calculating the image shifts. The subfield size was not changed. 



Removal of outliers 



Shift measurements that fail are considered as outliers. Sometimes when the 
shift calculations are made, the wrong set of granules is chosen for measuring 
the image offset. This means that the position in the reference subimagc is 
compared to another position in the other subimage. These positions do not 
correspond to the same granules and the shift measured will be completely 
wrong. Including these measurements when calculating the covariances strongly 
reduces the resulting Tq. 

A filter for outliers where therefore implemented in the code, before the covari- 
ances arc calculated. The standard deviation for all image shifts (for a particular 
data set) is calculated. Subfields with a larger shift than 3 standard deviations 
are considered as outliers and ignored in the covariance calculations. 
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th Chapter 

Results 



Data from all observing days in 2009 and for observing days to the 13th of 
September in 2010 have been run through the different inspection scripts cre- 
ated. Statistics presented in this chapter are made from these days. 

WFWFS data from ten days in 2010 have been fully processed. Five of these 
days (9th, 21st, 22nd and 28th of June as well as 16th of July) have been 
compared with data from the SHABAR. Comparisons, results and discussion 
from the first day, the 9th of June, are presented in this chapter. 



Inspection statistics 



Figure 6.1 shows how many image frames that were assigned with each status 
code. Figure [6?2| shows the same but for error status codes. The left histogram 
of each figure shows statistics from 2009 and the right histogram for 2010. The 
corresponding status codes can be found in Table |5. 1 | and a description of how 



to get the "combined" status codes are given in Sec. 5.1.1 



Status 2009 2010 

Good 853 171 (28 %) 292 581 (35 %) 
Warning 412 133 (13 %) 161 437 (19 %) 
Bad 1 792 330 (59 %) 377 365 (45 %) 
3 057 634 831 383 



Table 6.1: Number of files for each status good, warning and bad. 



A number of new (combined) status codes were introduced in 2010 because of 
the threshold change in status code 64. Some statuses were assigned in 2009 but 
haven't been assigned to any of the frames taken until the 13th of September 
2010. These are mostly statuses which are combinations of statuses assigned 
when the intensity is low. These statuses become more common in the end of 
the observation season. 
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Statuscodes for warning files 2009 



Statuscodes for warning files 2010 
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Figure 6.1: (a) Number of frames for each warning status code for data taken in 2009. 
(b) Number of frames for each warning status code for data taken in 2010. Status codes 
that occur in 2010 but didn't in 2009 are combinations with status 64 that didn't exist 
before due to the threshold change. Some statuses combined of statuses assigned to 
frames with low intensity didn't exist for frames in 2010 when the inspections where 
done. These might occur later in the observation season. 



statuscodes for "bad" files 2009 



Statuscodes for "bad" files 2010 
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Figure 6.2: (a) Number of frames for each error status code for data taken in 2009 
(b) Number of frames for each error status code for data taken in 2010. Status codes 
that occur in 2010 but didn't in 2009 are combinations with status 64 that didn't exist 
before due to the threshold change. 



As written in Sec. |5.1.3| the intensity of had frames from 2009 was checked to 
see how close to the threshold they were and when they were taken. Almost 
all frames with lower subaperture intensity than 1600 were taken in October or 
later. Intensity values between that and the threshold (2000) were spread over 
the summer. 
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Frames from 2009 that were deleted because of high RMS comes from two days, 
the 9th of July and the 24th of September. All other frames with high RMS 
are bad because of other things. Frames from the 9th of July have RMS higher 
than 20 but frames from the 24 of July have RMS values from 10 to 12, which 
is quite close to the threshold. These can not be examined by eye since all bad 
frames are deleted. 



^^0 



The correlation coefficients obtained when vq from different heights (measured 
the 9th of June) were correlated with each other, are presented in correlation 
matrices. Fig. 6.3 The correlation coefficient is described in Eq. (5.5) and 



Sec. 5.3 In Figure 6.3a the correlation matrix, for all values of tq, is shown. 

It can be seen that some correlation exists between layers at higher heights as 
well as between the lowest layer and these at higher heights. That the smallest 
values for tq at the highest heights were associated with poor ground layer seeing 
was found by Scharmer and van Werkhoven ppj . 

The correlation matrix changes when < 7.5 cm are excluded. The new corre- 



lation matrix is shown in Fig. 6.3b The correlation between heights (especially 
higher) are now less than in the previous case. The difference was first noticed 
when the limit was set to 7.5 cm. This limit was therefore used in the rest of 
this thesis. 





(a) (b) 

Figure 6.3: (a) Correlation matrix of ro for WFWFS layers. Black pixels indicate 
no correlation and white pixels show perfect correlation (when a layer is correlated 
with itself). The gray pixels show that there exist some correlation for layers at higher 
heights, (b) Correlation matrix of ro for WFWFS layers where the integrated ro < 7.5 
cm have been excluded. The correlation between higher layers are much less when 
these values are excluded. 



The calculated ro values, from the 9th of June, for the 7 different WFWFS 
heights are shown in Figs. |6.4f|6.6[ 

The high peaks for tq shown in Fig. |6.5b| are probably due to no turbulence 
region in this WFWFS layer. 
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WFWFS r„ at /i = m 
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Figure 6.4: ro calculated from WFWFS measurements at /i = m. 



WFWFS at )i=500 m 





08.00 10.00 



12.00 14.00 
Time 



(a) 500 m 



WFWFS at )i = 1300 m 



14. 



jLm 



08.00 10.00 12.00 14.00 16.00 

Time 



(b) 1300 m 
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Figure 6.5: ro calculated from WFWFS measurements at different heights 
h = [500, 1300, 3200, 5800] m. 
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WFWFS r„ at /i = 10800 m 



WFWFS r„ at /i =20000 m 
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Figure 6.6: ro calculated from WFWFS measurements at different heights 
h = [10800, 20000] m. 



WFWFS ro COMPARED WITH SHABAR ro 



The correlation matrix of tq for the 10 SHABAR layers first used (described 



in Sec. 5.3), Fig. 6.7a shows that there were too many SHABAR layers. The 



measured rg values in the different layers correlate strongly with each other, 
which indicates that the measurements are not independent. The layers were 
therefore changed so that the correlation between them were as low as pos- 
sible. The new boundaries between layers were chosen to be located at 1 
m, 10 m, 30 m, 200 m, 800 m, 3 km and 20 km. These 7 new heights de- 
fine boundaries between 6 new SHABAR layers which have average heights 
hshabar = [5, 20, 115, 500, 1900, 11500]. The corresponding correlation matrix is 
shown in Fig. |6.7b 



The 6 new SHABAR layers were compared with the 7 WFWFS layers. The 
correlation between rg calculated from WFWFS measurements at h — Q m, 
plotted against the best fitted tq from SHABAR data is shown in Fig. |6.8[ 
The correlation coefficient is 0.95. How much each of the SHABAR layers 
contributed to the best fitted tq for that WFWFS height {h ~ m) is shown 
in Fig. |6.9| This shows that the contribution from SHABAR data to the best 
fitted ro comes from layers below 200 meter. 



Figure 6.10 shows the correlation between WFWFS measurements and best 
fitted SHABAR rg for WFWFS-height h = 500 m. The correlation coefficient is 
0.27 which means that the correlation is not as strong as for the first layer but 
still there. The corresponding contribution from the different SHABAR layers 
is shown in Fig. |6.11[ which shows that the highest contributions come from the 
two layers located at h — 500 m and h = 1900 m. 
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(a) (b) 

Figure 6.7: (a) Correlation matrix of ro for tiie 10 SHABAR layers first used. Black 
pixels indicate no correlation and white pixels show perfect correlation (when a layer 
is correlated with itself). The different white-gray pixels show how much correlation 
that exist between layers, (b) Correlation matrix for the 6 SHABAR layers chosen to 
minimize the correlation between layers. 

The kernels for the SHABAR inversion converge at heights larger than 1 km 
[TTll25j which should set the limit for the SHABAR range. This means that the 
SHABAR most likely cannot distinguish layers at higher heights. 

No comparison between WFWFS height 1300 m and SHABAR data is shown 
in the figures below since these tq values were too high. No turbulence regions 
seem to be located at this height. 

It is hard to see any correlation between the WFWFS and the SHABAR for 
h — 3200 m, where the correlation coefficient is 0.07. The contribution to the 
fitted ro still comes mostly from expected layers {h = 500, 1900), even though 
the correlation is low. 

The correlation coefficient is even lower, -0.01, for h = 5800 m. The contri- 
bution from the lowest layer gets stronger for this height but the peak is still 
at higher heights [hghabar = 500,2000). Almost no correlation is found for 
WFWFS height h = 10800 m either, where the correlation coefficient is 0.16. 
The contribution to the best fit comes mostly from the first layer. 

Some correlation seems to be present again at higher WFWFS heights. Espe- 
cially for h = 20000 m, where the correlation coefficients is 0.22. This effect 
is however not real. The seeing is worse here than for the closest previous lay- 
ers and the small values of ro come from the poor ground layer seeing. The 
correlation found is because of the correlation at lower heights. 
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Figure 6.8: Correlation between best fit ro from SHABAR data and ro calculated 
from WFWFS measurements at ft = m. Correlation coefficient 0.95. 




Figure 6.9: Contribution to best fit ro for WFWFS height h = 0. 
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WFWFS r„ to fitted SHABAR r. 




0.2 0.4 0.6 0.8 1.0 1.2 1.4 

WFWFS r„ (with 7.5 cm r„ limit) for ft =500 [m] 



Figure 6.10: Correlation between best fit ro from SHABAR data and ro calculated 
from WFWFS measurements aX h = 500 m. Correlation coefficient 0.27. 
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Figure 6.11: Contribution to best fit ro for WFWFS height h = 500. 
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WFWFS r„ to fitted SHABAR r. 
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Figure 6.12: Correlation between best fit ro from SHABAR data and ro calculated 
from WFWFS measurements at ft = 3200 m. Correlation coefficient 0.07. 
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Figure 6.13: Contribution to best fit ro for WFWFS height h = 3200. 
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WFWFS r„ to fitted SHABAR r. 
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WFWFS r„ (with 7.5 cm r„ limit) for/i=5800 [m] 



Figure 6.14: Correlation between best fit ro from SHABAR data and ro calculated 
from WFWFS measurements aX h = 5800 m. Correlation coefficient -0.01. 
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Figure 6.15: Contribution to best fit ro for WFWFS height h = 5800. 
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Figure 6.16: Correlation between best fit ro from SHABAR data and ro calculated 
from WFWFS measurements at ft = 10800 m. Correlation coefficient 0.16. 




Figure 6.17: Contribution to best fit ro for WFWFS height h = 10800. 
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WFWFS r„ to fitted SHABAR r. 
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WFWFS r„ (with 7.5 cm r„ limit) for/i=20000 [m] 



Figure 6.18: Correlation between best fit ro from SHABAR data and ro calculated 
from WFWFS measurements aX h = 20000 m. Correlation coefficient 0.22. 
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Figure 6.19: Contribution to best fit ro for WFWFS height h = 20000. 
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Smaller subfield size and removed outliers 

The resulting covariances achieved after changing the circular subfield mask (ap- 
plied when calculating the image shift) are compared to the original covariances 
in Figs. [6:20p:2T| 

Changing the mask size from 16 to 12 pixels increased the number of failed 
shifts (larger than 3cr) from 1.35 % to 3.29 %. The percentage of failed shifts 
for the 8 pixel mask was 7.43 %. 

Longitudal covariances 

12| ^ ^ ^ ^ ^ ^ ^ ^ 1 

10- 




16x16 



Longitudal covariances, outiiers removed 




16x16 
(a) 

Figure 6.20: Resulting covariances obtained with a 16 x 16 circular subfield mask for 
image shift calculations compared with the resulting covariances when a 12 x 12 circular 
subfield mask was applied. The different colours show different angular separation, 
(a) Covariances calculated with all shift measurements. The blue points that have a 
much higher covariance when using the 12 pixel mask corresponds to no subaperture 
separation s — 0. (b) Covariances calculated after shift measurements larger than 3a 
was removed. 
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Both the 12 and 8 pixel mask showed lower covariances than the 16 pixel mask, 
before the outliers where removed. This was stiU the case for the 8 pixel mask 



after the outliers where removed, see Fig. 6.21b As can be seen in Fig. 6.20a 



the 12 pixel mask and the 16 pixel mask correlate well after removal of the failed 
shifts. 



Longitudal covariances 




(a) 



Longitudal covariances, outliers removed 




(b) 

Figure 6.21: Resulting covariances obtained with a 16 x 16 circular subfield mask for 
image shift calculations compared with the resulting covariances when a 8 x 8 circular 
subfield mask was applied. The different colours show different angular separation. 

(a) Covariances calculated with all shift measurements. The blue points that have a 
much higher covariance when using the 12 pixel mask corresponds to no subaperture 
separation s = 0. (a) Covariances calculated with all shift measurements. 

(b) Covariances calculated after shift measurements larger than 3a was removed. 



48 




th Chapter 



Conclusions 



Some conclusions drawn from the results in Chapter [6] are presented in this 
chapter together with some thoughts about next step and future improvements. 



ro AT DIFFERENT HEIGHTS 

Scharmer and van Werkhoven [20. found that the smallest values for tq at the 
highest height were associated with poor ground layer seeing. The images get 
noisy and the cross-correlations fail if the small values of ro is kept. They 
determined the value for when rg should be rejected to 7.5 cm, which is about 
75 % of the subaperture size. 

As discussed in Sec. |6.2| and shown in Fig. |6.3| the same result was found during 
this project. It was therefore concluded that this limit shall be used. 



The values of ro differs a lot from height to height, which is shown in Figs. 6.4 



6.6 ro is small for the first height, ft, = m, as well as for the highest height 
h = 20 km. It can therefore be concluded that one turbulent layer is located at 
the ground level and one at higher altitude. 

The intermediate layers have rg larger than 40 cm and ro is sometimes really 
large for h = 1300 m. It can be concluded that the turbulent layer at this height 
is weak. 



Correlation between WFWFS and SHABAR 

The correlation between calculated ro from the WFWFS and from the SHABAR 
is best for lower layers, h = and h = 500 m. Next WFWFS node h = 1300 
showed very varying and high ro and no comparison with the SHABAR was 
therefore possible. 

The small correlation found for higher layers still have most contribution from 
layer hghabar — 500 m. It therefore seems likely that the SHABAR cannot 
distinguish layers at these heights. This is consistent with the convergence of 
the SHABAR kernels 
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CHAPTER 7. CONCLUSIONS 



New subfield mask and outliers removed 



The comparison between the covariances achieved with different sizes of the 



circular subfield mask applied when calculating the image shifts, Figs. 6.20 



6.21 showed that the 12 pixel mask is useful but the 8 pixel mask is not. 



Removal of outliers is an improvement that will be used in the further analysis. 



Future 

Even though data from several days are reduced it is only one day (2010-06-09) 
that is compared in detail with the SHABAR and this comparison was done 
before the outliers were removed from the WFWFS data. A new comparison, 
without these outliers, would be interesting. Comparing more days both in the 
beginning and end of the observation season would also be interesting. The 
reduction script can, after being started, work autonomously with reduction of 
data from a whole year. The problem is where to store the reduced data. The 
resulting covariance map takes up little space but the steps in between produce 
quite large files. One days reduction takes up 20-170 GB depending on the 
number of good sets and the number of references used in the reduction. 

A similar WFWFS will be placed at the Vacuum Tower Telescope (VTT) on 
Tenerife. After comparison of the data from the two WFWFS, the different sites 
can be evaluated in terms of strength, location and stratification of the seeing. 

Comparison with a long baseline SHABAR and a WFWFS (at both sites) will 
also be done. This SHABAR will have a baseline of 3.5 m. 

After all these further reductions and comparisons, it might be possible to decide 
where to build the European Solar Telescope and how to make a good design 
for its Multi Conjugate Adaptive Optic system. 
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8 



th Chapter 



Author's 



CONTRIBUTIONS 



This chapter is a short summary of what I have done during my master project. 
Inspection 

I changed some thresholds for the inspection to easier distinguish between bad 
and good frames. I also wrote several inspection scripts that makes it easier to 
compare the results from the automated inspection. 

All data from 2009 and most of 2010 were inspected using these scripts and the 
result can now be used to make a decision of which files to save. 



Reduction in one go 

I wrote a reduction script so that all the reduction steps could be evaluated 
after each other. It is now possible to reduce data from a whole year in one go. 

Change of subfield mask 

I changed the circular subfield mask used when calculating the image shifts to 
test if better lieiglit resolution could be achieved. Two new sizes where tested. 



Removing outliers 

The S-DIMM+ analysis was changed so that outliers (failed shift measurements) 
are removed before the covariances are calculated. 



Comparison between the WFWFS and the SHABAR 

I compared the results from the Wide Field Wavefront Sensor with the Shadow 
Band Ranger mounted on the SST. 
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CHAPTER 8. AUTHOR'S CONTRIBUTIONS 



Documentation 

I have documented the changes and updated the existing documentation for 
Astooki written by T.I.M. van Werkhoven [27] , 

Updating software 

I have had discussions with T.I.M. van Werkhoven about implementing my 
changes into his software permanently. 
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Inspection scripts 



All inspection scripts store the output to a text file after every step. Not needed 
if you're only interested in the output of the last step but it turned out to be 
quite handy to be able to look at all of them at different times. 

Shown below is a small bash script that inspects which status codes all bad 
frames had. The path to all log files from 2010 are stored in logfiles.txt. 

Code snippet A.l: Inspection of bad frames (inspection. sh) 

# 

# Locate logfiles and sort out files with only bad data 

# 

# Store in a badlog 

logf iles = ' cat logfiles.txt' 
for logfile in $logfiles 
do 

bad='cat $logfile I grep -c BAD' 

if [ $bad != ] 

then 

echo $logf ile : $bad >> badlog.txt 

f i 
done 

# 

# Read path to logfile from badlog.txt 
FILES='awk -F '{print $!>' badlog.txt' 

# 

# Find status codes and store in a file 

# 

for FILE in $FILES 

do 

# Display status for files with error status 

grep BAD $FILE I cut -d',' -f 4 I sort -u I cut -d ' - ' -f 2 

done>statuses .txt 

# Make a sorted list of available status codes 
sort -nu statuses.txt >statuscodes_bad . txt 
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APPENDIX A. INSPECTION SCRIPTS 



rm statuses.txt 



# 

# List all directories (data sets) and number of bad files 

# 



# Loading the statuses into an "array" 
STATUSCODES= ' awk '{print $1}' statuscodes_bad.txt' 

for FILE in $FILES 
do 

# Display which directory your looking in 
id='echo $FILE I cut -f 4 -d/ ' 

echo "Processingudirectory:u$id" 

# Display status for files with bad status 

grep BAD $FILE I cut -f 4 -d',' I cut -d'-' -f 2 > 
filestatuses . txt 



for code in $STATUSCODES 
do 

count='grep -wc "$code" filestatuses.txt' 
echo " Number u of uf iles ywithuStatusuScode : □$ count " 
done 

done >bad . txt 

rm filestatuses.txt 

# 

# Read the bad files and store them in an array 

# 

# Number of statuscodes 

columns=$(( 'wc -1 statuscodes_bad.txt I cut -f 1 -d' ' ' + 
1 )) 

rows=$(( 'wc -1 badlog.txt I cut -f 1 -d ' ''+!)) 

row = l 
column=l 

while [ $row -le $rows ] 
do 

while [ $column -le $(( $row * $columns )) ] 
do 

text='head -n $column bad. txt I tail -n 1 I cut -f 2 -d 

' : ' I cut -f 2 -d ' ' ' 
printf "Z.SuXt" $text 
let column++ 
done 
echo "" 
let row++ 
done >badarray . txt 
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Reduction script 



The code snippet in |B.l| is a part of the bash script reduce. sh, that calls and 
sends input variables to pyatk.py. pyatk.py is the base class for all other 
"tools" provided by Astooki. 

The different steps in the reduction described in Sec. |5.2| are shown below. The 
beginning of the script, where flat fields and dark frames are chosen, are not 
shown below. Neither are the loops that make the reduction run autonomously 
for all specified days. 

Code snippet B.l: Reduction script (reduce. sh) 



# Subimage mask 

# Equal for all data sets (runs) 

pyatk.py samask -vv -d / scrat ch/ wf wf s /proc / simask \ 
--file simask. CSV \ 
--rad 1024 \ 
--shape circular \ 
--sasize 175,154 \ 
— pitch 194,166 \ 
--xoff 0,0.5 \ 
--disp 1017,1000 \ 
--scale=l \ 
--plot 



# Optimizing the subimage mask for that set 

pyatk.py saopt -vv -d / s cr at ch / wfwf s /proc / $day/ $imset - simask 
-fopt \ 

--mf / s cr at ch/ wf wf s /pr oc / s imask/ s imask - or igin . cs V \ 
— ff /scratch/wfwf s/$f f day/$f f day-$f f set/$f fbest — fm 300 
\ 

--file simask - or ig - f opt . CSV \ 
-- saif ac 0.8 \ 
--rad 1024 \ 
--plot 
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APPENDIX B. REDUCTION SCRIPT 



# Get sasize values (from prev . step) 

sasize='grep subimage / scr at ch/ wf wf s /proc / $day / $imset - simask 
-f opt /astooki -log I cut -f 3 -d '(' I cut -f 1 -d 
sort -u ' 

# Subfield mask 

pyatk.py sfmask -vv -d / s cr at ch/ wf wf s /pr oc / $day / $imset - 
sfmask \ 
--file sfmask -big . CSV \ 
--sf size=113 ,89 \ 
-- sas ize= $sas ize \ 
-- over lap =0 , \ 
--border=30 ,30 \ 
--plot 

# 

# Static offset 

pyatk.py shifts -vv -d / s cr at ch/ wf wf s /pr oc / $day / $imset - 
statoff \ 

--ff /scratch/wfwf s/$f f day/$f fday-$f f set/$f fbest — fm 300 
\ 

--df /scratch/wfwf s/$ddday/$ddday-$ddset/$ddbest --dm 300 
\ 

--safile /scratch/wfwfs/proc/$day/$imset-simask-fopt/ 

simask -orig - f opt . CSV \ 
--sffile /scratch/wfwfs/proc/$day/$imset-sfmask/sfmask-big 

. CSV \ 

--range 7 --nref 1 --mask none --plot \ 
/scratch/wfwfs/$day/$day-$imset/wfwfs_survey_im* 

# 

# Updated subimage mask (with offset) 

pyatk.py saupd -vv -d / scr at ch/ wf wf s /proc / $day / $imset - 
statoff \ 

--mf /scratch/wfwfs/proc/$day/$imset-simask-fopt/simask- 

or ig - f opt . CSV \ 
--offset /scratch/wfwfs/proc/$day/$imset-statoff/static- 

offsets.csv \ 
--plot 

# 

# Subfield mask with smaller subfields 

# Pick the smallest sasize (from prev. step) 

xsasize=' cat /scratch/wfwfs/proc/$day/$imset-statoff/simask- 
or ig - f opt -updat ed . cs V I head -n 3 I tail -n 1 I cut -f 1 
-d 

ysasize='cat /scratch/wfwfs/proc/$day/$imset-statoff/simask- 
or ig -f opt -updat ed . CSV I head -n 3 I tail -n 1 I cut -f 2 
-d , ' 

pyatk.py sfmask -vv -d / s cr at ch/ wf wf s /pr oc / $day / $imset - 
sfmask \ 
--file sfmask -16x16 . CSV \ 
-- sf size =16 ,16 \ 
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64 --sasize=$xsasize,$ysasize \ 

65 --overlap =0.7, 0.7 \ 

66 --border = 7,7 

67 

68 # 

69 # Measure the subfield shifts 

70 pyatk.py shifts -vv -d / scr at ch/ wf wf s /proc / $day / $imset - 

subshift -16x16 \ 

71 — ff /scratch/wfwf s/$f f day/$f f day-$f f set/$f fbest — fm 300 

\ 

72 --df /scratch/wfwf s/$ddday/$ddday-$ddset/$ddbest --dm 300 

\ 

73 --safile / s cr at ch / wf wf s /pr oc / $day / $imset - st at of f / simask - 

or ig -f opt -updated . CSV \ 

74 --sffile / s cr at ch / wfwf s /proc / $day/ $imset - sf mask/ sf mask - 16 

xl6 . CSV \ 

75 --range 7 --nref 2 --mask circular \ 

76 / scratch/wfwfs/$day/$day-$imset/wfwfs_survey_im* 

77 

78 # 

79 # Subaperture mask for the lenslets 

80 # Equal for all runs 

81 pyatk.py samask -vv -d / scr at ch/ wfwf s /proc / samask \ 

82 --file samask-ll.csv \ 

83 — rad 0.52 \ 

84 --shape circular \ 

85 — sasize 0.098,0.098 \ 

86 — pitch 0.098,0.0849 \ 

87 — xoff 0,0.5 \ 

88 --disp 0,0 \ 

89 --scale = l --plot 

90 

91 # 

92 # S-DIMM+ analysis 

93 pyatk.py sdimm -vv -d / s cr at ch / wfwf s /proc / $day/ $imset - sdimm 

-16x16 \ 

94 --skipsa -1 \ 

95 --shifts / s cr at ch / wfwf s /proc / $day/ $imset - subshift - 16 xl6 / 

image - shift s . npy \ 

96 --safile / s cr at ch / wfwf s /proc / samask/ samask -11 - c entro id . CSV 

\ 

97 --sffile / s cr at ch / wfwf s /proc / sf mask/ sf mask - 16 xl6 . CSV 
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APPENDIX B. REDUCTION SCRIPT 
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